/** 
 * \file gddot.cu
 * \author Kyle E. Niemeyer
 * \date 10/04/2011
 *
 * Based on "ddot.f" from BLAS.
 *
 */

////////////////////////////////////////////////////////////////////////

/** gddot forms the dot product of two vectors.
 * Uses unrolled loops for increments equal to one.
 *
 * \param[in] n       vector dimension
 * \param[in] dx      array, dimension n
 * \param[in] incx    increment between elements of dx
 * \param[in] dy      array, dimension (n-1)*|incy|+1
 * \param[in] incy    increment between elements of dy
 * \return    ddot    dot product
 */
__device__ __inline__ double gddot ( int n, const double *dx, int incx, const double *sy, int incy )
{
  
  double ddot = 0.0;
  
  if ( n <= 0 ) return;
  
  if ( ( incx == 1 ) && ( incy == 1 ) ) {
    
    // both increments equal to 1
    
    // clean-up loop
    
    uint m = n % 5;
    
    if ( m != 0 ) {
      for ( uint i = 0; i < m; ++i ) {
        ddot += dx[i] * dy[i];
      }
      
      if ( n < 5 ) return ddot;
    }
      
    for ( uint i = m; i < n; i += 5 ) {
      ddot += dx[i] * dy[i]
            + dx[i + 1] * dy[i + 1]
            + dx[i + 2] * dy[i + 2]
            + dx[i + 3] * dy[i + 3]
            + dx[i + 4] * dy[i + 4];
    }
      
  } else {
    
    // increments not equal to 1 or unequal increments
    
    int ix = 1;
    int iy = 1;
    
    if ( incx < 0 ) ix = (-n + 1) * incx + 1;
    if ( incy < 0 ) iy = (-n + 1) * incy + 1;
    
    for ( uint i = 0; i < n; ++i ) {
      ddot += dx[ix] * dy[iy];
      ix += incx;
      iy += incy;
    }
    
  }
  
  return ddot;
  
}